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Abstract. Many challenging tasks in sensor networks, including sensor 
calibration, ranking of nodes, monitoring, event region detection, collab- 
orative filtering, collaborative signal processing, etc., can be formulated 
as a problem of solving a linear system of equations. Several recent works 
propose different distributed algorithms for solving these problems, usu- 
ally by using linear iterative numerical methods. 

In this work, we extend the settings of the above approaches, by adding 
another dimension to the problem. Specifically, we are interested in self- 
stabilizing algorithms, that continuously run and converge to a solution 
from any initial state. This aspect of the problem is highly important 
due to the dynamic nature of the network and the frequent changes in 
the measured environment. 

In this paper, we link together algorithms from two different domains. 
On the one hand, we use the rich linear algebra literature of linear it- 
erative methods for solving systems of linear equations, which are natu- 
rally distributed with rapid convergence properties. On the other hand, 
we are interested in self-stabilizing algorithms, where the input to the 
computation is constantly changing, and we would like the algorithms to 
converge from any initial state. We propose a simple novel method called 
SS-Iterative as a self-stabilizing variant of the linear iterative meth- 
ods. We prove that under mild conditions the self-stabilizing algorithm 
converges to a desired result. We further extend these results to handle 
the asynchronous case. 

As a case study, we discuss the sensor calibration problem and provide 
simulation results to support the applicability of our approach. 

1 Introduction 

Many challenging tasks in sensor networks, for example distributed rank- 
ing algorithms of nodes and data items [3], collaborative filtering [1], local- 
ization [10], collaborative signal processing [12], region detection [9], etc., 

* Part of the work was done while the author visited Cornell university. The work was 
funded in part by ISF. 
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can be formulated as a problem of solving a linear system of equations. 
Several recent works [10], [12], [9] propose different distributed algorithms 
for solving these problems, usually by linear iterative numerical methods. 

In this work, we extend the settings of the above approaches by adding 
another dimension to the problem. Specifically, we are interested in self- 
stabilizing algorithms, that continuously run and converge to a solution 
from any initial state. This aspect of the problem is highly important due 
to the dynamic nature of the network and the frequent changes in the 
measured environment. 

As a case study, we show that the calibration of local sensors' readings 
can be formulated as a linear system of equations Ax = b, where x 
represents the calibrated output reading, b represents the local reading, 
and A represents a weighted communication graph. However, our work 
is general and can be applied to any problem that can be formulated 
as a distributed solution to a linear system of equations, including the 
previous works mentioned above. 

Consider a distributed system of sensors measuring real-world data. 
Sensors are located in different areas; for example, the senors are spread 
throughout a building and they measure the temperature to adjust the 
heating or cooling. We would like the collected data to be as reliable as 
possible, reflecting closely the changing environmental conditions. One of 
the obstacles we face when designing algorithms that collect data from a 
sensor network are measurement errors. There are two main types of inac- 
curacies of sensors' measurements: noisy environment and sensing equip- 
ment which is not calibrated. It is desirable that sensors could execute 
a distributed algorithm for calibrating their environmental readings. In 
this setting sensors are allowed to communicate among themselves, using 
data from other nodes to affect their reported individual reading. Fur- 
thermore, we would like our calibration algorithm to have fault-tolerance 
properties. Specifically, we are interested in self-stabilizing algorithms [7] 
which converge to an optimal solution from any initial state. Observe that 
self-stabilization helps also in deploying the sensors. There is no need to 
explicitly synchronize the sensors, once enough of them are deployed and 
begin functioning the results will converge to the expected value. 

The main challenge we have faced in this work, is that in the classical 
linear algebra literature, b is assumed to be constant. In our settings, the 
environment is constantly changing and the computed algorithm never 
terminates, leading to constantly changing values of b. In this paper, 
we ask the following question: "Is it possible to devise a self-stabilizing 
numerical iterative method?" We answer affirmatively, and show that 
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under minor conditions it is possible to devise a self-stabilizing algorithm 
that solves a dynamic system of linear equations, where the input to the 
system is constantly changing. 

To the best of our knowledge, this is the first work tackling this chal- 
lenging problem. We believe that our approach can have numerous appli- 
cations in the field of distributed self-stabilizing computation. 

Other works discuss fault tolerance aspects of distributed computa- 
tion. For example, overcoming faults in sensors by averaging the input 
was investigated in [11] providing a centralized algorithm. Quantifying 
faulty nodes' effect on the system's output is discussed in [8] and [5]. 
These papers consider bounded input paths and their effect on the sta- 
bility of the output. In [6] infinite input paths are considered under the 
assumption that only specific sensors' input may change. All three pa- 
pers consider discrete input values, as opposed to a continuous set of 
input values discussed in this paper. 

The paper is constructed as follows. Section 2 defines the model and 
problem definition. Section 3 presents our novel algorithm SS-Iterative. 
Section 4 analyzes our algorithm and gives bounds on the convergence 
rate. Section 5 presents experimental results of running SS-Iterative 
using sample topologies. Section 6 extends our construction to the asyn- 
chronous case. We conclude in Section 7. 

2 Model and Problem Definition 

We model the sensor calibration problem as follows. Given a directed 
communication graph G = {V,E), V is the set of nodes V = {pi, . . . ,Pn}, 
E is the set of weighted edges connecting them (weights can be negative) 
and N{pi) denotes pi's neighbors. Edge weights are used to model the 
directional dependence between nodes' outputs; i.e., if Wp^^p. = then 
there is no edge between p, and pj and their output is not directly de- 
pendent on each other. In addition, we require a non-zero self connected 
edge, Wp^^p■ 7^ 0, which represents the weight of pi's own input. 

Initially, we assume a synchronous system: during a single round of 
communication, any pair of connected nodes may send a single message 
on each directed edge. Each round r, each node pi has a scalar input value 
Ipi{r), which represents the local reading of the sensor.'^ In addition, pi 
outputs its output value, which is denoted by Op^{r); both inputs and 
outputs are from the domain of real numbers. Denote by I{r) the input 

^ For simplicity of notations we use scalar variables in the paper. An extension to the 
vector case (where each sensor measures a set of measurements) is immediate. 



4 



vector of the entire system at round r, and by 0(r) the output vector of 
the system at the end of round r. In Section 6 we relax the assumption of 
synchronous rounds and provide a variant of the algorithm which works 
in asynchronous settings. 

The schematic operation of each node pi at round r is composed of 
the following steps: (a) read the value of Ipi{r); (b) send messages; (c) 
receive messages; and (d) do some processing and output Op. (r). Then a 
new round is started, and the nodes continue so forever. 

Definition 1. A configuration C of the system at round r consists of 
the state of each node prior to performing any operation at round r; this 
configuration is denoted by C{r). 

Definition 2. An input sequence Z of length i is a list of £ vectors 
such that each v ^ I is a possible input vector of the system (i.e., v € 2), 
the domain of allowed values). An output sequence O of length i is 
a list such that each v € O can potentially be an output vector of the 
system. 

Definition 3. A step from configuration C to configuration C on input 
vector V is legal if C is reached from C by the system when having v as 
the input vector, u is produced by a legal step if u is the output vector 
of the system resulting from such a legal step. 

Definition 4. A run of a system on input sequence I = {v(l), . . . , 
starting from configuration C{r) is the sequence C{r),0{r),C{r + l),0{r + 
1), . . . s.t. for any i > 0; the step from C{r + i) to C{r + i + 1) on input 
vector v(i + 1) is legal, and 0{r + i) is produced by that legal step. The 
system is said to produce the output sequence O = {0(r), . . . , 0(r+£— 1)}. 

In the special case when the sensor observations (the input to the 
system) are fixed, the output decision of the sensors should converge 
to a solution that preserves the linear relations among node inputs and 
outputs. More formally, consider an input sequence I of identical input 
vectors; i.e., Z = {v,v,v, ...}. It is desired that for such an input a 
run from any configuration C on I would end up producing an output 
sequence O = {u(l), u(2), . . . } such that ||u(i) — u|| ^ as z ^ oo, for 
a u that solves the following linear system of equations: 



P3&N{pi) 



(1) 



5 



We assume that the above equations are uniquely solvable, denoting u as 
the solution to v. 

One of the most efficient distributed approaches for solving a set of 
linear equations of the type Ax = b is by using linear iterative algo- 
rithms. Unlike Gaussian elimination, which has a cost of O(n^), where n 
is the number of variables, an iterative algorithm usually solves a system 
of linear equations in time of 0{v?r, ) where r is the number of itera- 
tions, which is typically logarithmic in n. These algorithms are naturally 
distributed and work well in asynchronous settings. Furthermore, when 
converging, the algorithms converge to a solution from any initial state. 
An excellent survey of such methods is found in [2]. 

The main novel contribution of this paper is in analyzing the self- 
stabilizing properties of algorithms from the linear iterative methods do- 
main. In a practical setting, it is highly unreasonable to assume that 
sensor readings do not change over time. However, it is reasonable to 
assume that at steady state the change in sensor readings is bounded. In- 
formally, in this work we show that once the input readings are bounded, 
the output solution is bounded as well. This useful observation enables 
us to tie together numerical iterative methods and dynamically changing 
environments in a self-stabilizing manner. 

The following definition bounds the change in sensor observations: 

Definition 5. An input sequence I = {v(l), v(2), . . . ,v(^)} is (5-bounded 

around v if for every i, 1 < i < i, it holds that ||v(i) — v||^ < 5.^ 

Definition 5 states that a sequence I is (^-bounded if all the vectors 
in I are bounded within an n dimensional hypercube with an edge 25, 
centered around a point v. We note that once changes in the input are 
not bounded, then no efficient algorithm (especially in a network that is 
sparsely connected) can calculate the output fast enough. For example, 
if the diameter of the communication graph is T>, for some system of 
equations it would take at least T> rounds for the information exchange 
for input readings at one side of the network to propagate to the other 
side of the network. 

Definition 6. Let X he an input sequence that is 5-bounded around v 
and let u be the solution to input v. A run from configuration C on input 
sequence I e-converges to its solution if the produced output sequence 
O = {u(l),u(2),... u(Z\t)} satisfies that | |u(Z\t) - u)| |^ < e(Z\t,5,C); 
where e is a function of At, 5 and C. 

l|x|loo = maxi{|x,|}. 
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Definition 6 requires that if - starting from configuration C - the inputs 
are in an n dimensional hypercube of radius 5 around v then the output 
at time At is bounded within some n dimensional hypercube around u 
with radius €{At,d,C). We aim at an e{At,5,C) that decreases as At 
increases, as long as the inputs are bounded by the same v-centered, 6- 
radius hypercube. Clearly, for 6 > 0, e{At, 5,C) > for any At. That is, 
there is some minimal radius 6' > around u s.t. we cannot ensure a 
tighter bound. 

The above definition considers a single initial configuration, and a 
single input sequence I. We are interested in an algorithm that works for 
all initial configurations and all input sequences. 

Definition 7. An algorithm A e-converges for 5 -bounded input sequence 
I if every run (from any configuration) on I, e-converges to its solution. 
An algorithm A e-always converges if for every 6-bounded input se- 
quence I, A e-converges. 

Definition 7 formally defines the problem at hand, as an algorithm 
A that always converges has the desired self-stabilizing property: for any 
system state, once the sensors' readings changes are bounded, the change 
in output of the entire system is bounded as well. 

Our goal is to find an algorithm A that is e-always converging for 
a provably "good" e. Moreover, we aim at having A efficient also in its 
message complexity and simplicity of code, allowing lightweight sensors 
to actually implement it. 

3 Our Proposed Solution 

An equivalent formulation of the update rule Eq. (1) is 



The above equation states a condition on pj's output, in regard to pi's 
inputs and pi's neighbors' output. Thus, it encapsulates the requirement 
that different nodes influence each other's reported readings, while taking 
into consideration their local readings as well. 

Since Wp^^p^ ^ 0, the above equation can be stated as: 
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By denoting Wij = -^/^ (for i j) and Wi^i = we get: 

j 

Let be the matrix that has Wij as entries, Eq. (2) can be written 
in hnear algebra notation, (s.t. it appUes to all nodes simultaneously): 

Wu = v . (3) 

If we consider a non-self-stabilizing system in which the inputs do 
not change (that is, the input is fixed to v), then Eq. (3) can be seen as 

= b, where A and b are given. In such a case, we are interested in 
finding the value of x, which is a vector of n unknown variables. However, 
we are interested in the case where v changes over time, and thus Eq. (3) 
does not describe the problem properly, but rather helps in understanding 
the motivation for our solution. 

We use a modified update rule (relative to Eq. (1)): 

Op, (r + 1) = Wp^,p^ ■ Ip^ {r + l) + J2 Wp^,P, ■ Op, (r) . (4) 

Clearly, for the case of 5 = 0, a 0-bounded input sequence X, if {Op^ {r+ 
1) — Op^{r)) — > as r ^ oo then Eq. (4) converges to the solution of 
Eq. (1). Thus, if the update rule of Eq. (4) is executed simultaneously 
by all nodes, and for all of the nodes {Op-{r -|- 1) — Op^(r)) — > 0, then 
it also solves Eq. (3). That is, if each node locally executes Eq. (4) then 
the global solution is reached. This observation motivates algorithm SS- 
Iterative in Figure 1. 

Remark 1. In SS-lTERATlVE there is no notion of the "current round 
number r". That is, pi reads and writes to the variables Ip. and Op. 
without being "aware" of r. When we discuss the algorithm "from the 
outside", we will consider /pj(r) and Op^{r) instead of just /p.. Op.. 

Consider pi is running at round r + 1. When pi performs Line 03, it 
sends the value of Op.. The last time Op. was updated was at Line 04 and 
Line 06 of round r. Thus, the value sent by pi at round r -|- 1 is actually 
Op^{r). Therefore, the values received from pj by pi and used to update 
Op. {r + \) are Op. (r) . However, the value read by pi in Line 04 is the value 
of ipi(r + 1). Concluding that pi updates Op^(r -|- 1) exactly according to 
Eq. (4). 

Remark 2. Each node pi must know the values of Wp^^p■ as "part of the 
code". Thus, these values cannot be subject to transient faults. 
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Algorithm SS-Iterative 

01 : Each round do: /* executed on node pi */ 

/* send current value of Op- to all neighbors */ 
02: for each pj G N{pi) 
03 : send Op- to pj ; 

/* update Op- according to values sent by neighbors */ 
04: set Op. := Wp.,p. ■ Ip.\ 
05 : for each value Op^ received: 
06: update Op. := Op. + Wp^^p- ■ Op - ; 

07: od. 

Fig. 1. A self-stabilizing iterative algorithm. 



4 Analysis of SS-Iterative 

[2] shows that the update rule Eq. (4) can be written in linear algebra 
form as 

0{r + l)=AI{r + l) + BO{r), (5) 

where A is a diagonal matrix with Wp.^p. in the main diagonal, and Bij = 
Wp^^p. for i / j 

A^{diag{W})-^ , B^-{AW-ln^n) (6) 

where Inxn is the identity matrix. Using this update rule to solve a set 
of linear equations iteratively is known as the Jacobi algorithm. 

As noted in Section 2, when the input sequence is constant {i.e., 
I{r) = V for all r) the iterative execution of the above equations converges 
to u = + i?u, which is the same as u = VF~^v, thus solving Eq. (3). 
Following, we analyze the result of iteratively applying these equations 
for (5-bounded input sequences. 

Let X be an input sequence of length i that is (5-bounded around 
vector V. That is, X = I{r),I{r + 1), . . . , /(r + £ — 1) for some round 
r. Note that SS-Iterative saves a single scalar variable at each node, 
and thus the configuration of round r + 1 can be defined by the value of 
0(r) at round r. Consider SS-Iterative's run, starting from an arbitrary 
configuration at round r. We aim at showing that 0{r-\-At) is bounded by 
a hypercube centered at u. Denote by c{At) = 0(r -|- At) — u. If we show 
that ||c(Z\t)||j^ is bounded (as At increases), then 0{r + At) is within a 
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bounded hypercube centered at u. Consider c(l): 

c(l) = 0(r + l)-u 

= AI{r + 1) + BO{r) - {Av + Bu) 

= A{I{r + 1) - v) + B{0{r) - u) 

= A{I{r + 1) - v) + Bc{0) . (7) 

Since T is a 5-bounded input sequence around v, each I{r + At) can 
be denoted as v + P(r + At) s.t. V{r + At) G M" is a vector, and 
\\V{r + At)\\^ < 6. That is, V{r + At) = I{r + At) - v. 

Claim. At round r + At, it holds that c{At) = Efio ^ B^AVir + At - 
j) + S^*c(0). 

Proof. Proof by induction. The base of the induction was shown for c(l); 
see Eq. (7). Assume that the claim holds for At = k. Thus, c{k) = 
Ej=d B^AV{r + k-j) + B''c{0). By the update rule in Eq. (5), we have 
that 0(r + A; + 1) = AI{r + A; + 1) + BO{r + k). Combining the two 
equations implies 

c{k + 1) = 0{r + A; + 1) - u 

= AI{r + k + l)+ BO{r + k) - {Av + Bu) 
= A{I{r + k+l)-v) + B{0{r + k) - u) 
= AV{r + k + 1) + Bc{k) 

k-l 

= AV{r + A; + 1) + ^ 5-'+^y4P(r + k - j) + 5'=+^c(0) 

j=0 
k 

= AV{r + A; + 1) + ^ B^AV{r + k + 1 - j) + S'=+^c(0) 

k 

= B^AV{r + k + l-j)+ S*^+^c(0) . 
i=o 

Thus, if the claim holds for At = A; it also holds for At = k + 1; and we 
have that the claim holds for all At > 0. □ 

Definition 8. A matrix Mnxn is diagonally dominant if\Mii\ > Ej-_^-\Mij 

A matrix Mnxn is normalized diagonally dominant (normalized, for 
short) if M is diagonally dominant, and |Mjj| > 1. 
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Lemma 1. For a normalized diagonally dominant matrix W , it holds 

that \\A\\^ < 1 and ||-B||^ < 1, where A,B are defined in Eq. (6) and 
^ . 1 1 /I „i I 



1 



Proof. A is zero except for its main diagonal for which Ai^i 
Since \Wii\ > 1, we have that |A 
||x||^. Furthermore, maxx^o Trrr 

II M 

— '^Pi,Pj for i j and for i — j. Since W is assumed to be 
normaUzed diagonally dominant, we have that X^j^j \ < thus 
Y,jj^i\wpi,Pj\ < 1- Therefore, Y,j\Bij\ = < 1 all i. In 

total, for any x we have ||i3x|| < ||x|| , leading to ||5||„ < 1- □ 



t,i\ ^ 1- Thus, it holds that ||Ax||^ < 
< 1, i.e., Halloo < 1- Regarding B, 



If is a diagonally dominant matrix then node pj's own input effects 
Pi's output more than the sum of all of pj's neighbors outputs. That is, 
the weight of pi's input is at least the sum of weights of p^'s neighbors 
outputs. 

Theorem 1. Given a normalized diagonally dominant and invertihle W , 
there are constants ci,C2, where ci > 0, and 1 > C2 > 0, such that SS- 
Iterative e-always converges with e(Z\t, 5, C) = 5-ci+{c2)^*-\\0{r) — u||^ 

Proof. By Lemma 1 it holds that < 1 and < 1- Consider 

a (5-bounded input sequence I around v, and SS-Iterative's run start- 
ing from an arbitrary state 0(r). We are interested in the behavior of 

l|c(^t)lloo: 



B^AV{r + At- j) + B^*c(0) 

j=0 



< 



At~l 



B^AV{r + At 



j=0 
At-l 



+ ||S^*c(0) 



< Y \\B\U\AVir + At-j)\\^ + \\B\C'\\cm\, 



j=0 



<6-\\A\\^Yl II^IPoo + ll^ll^*l|c(0)||, 



j=0 



,At 



\-\B\r + 11-^11^' ii^wiic 



(8) 
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For an input sequence X that is 5-bounded around v, denote by u the 
solution to the original system of equations Ww = v. By Eq. (8), 

1 - II Rll^* 

\\0{r + At)-n\\^<6-\\A\\^-^^ + \\B\C\\cm\^. 

1 - Halloo 

Since < 1, we have that ^_||^|| < ^_||^|| — and by setting ci = 

^HjlHH' it holds that \\A\\^ iHlBlf - ^i- setting C2 = \\B\\^ and 
recalling that c(0) = 0(r) — u we are done. □ 

Theorem 1 states sufficient conditions s.t. SS-Iterative e-always con- 
verges. Moreover, the algorithm SS-Iterative is lightweight, as it re- 
quires nodes to send only a single value to every neighbor on each round. 



5 Experimental Results 

For illustrating the behavior of our proposed algorithm, we have simu- 
lated SS-Iterative using two sample topologies of one hundred nodes. 
Figure 2 depicts a circular topology where each node is connected to its 
left and right neighbors. Figure 3 shows a random unit disc graph, where 
nodes are randomly spread on a plane, and each node is connected to the 
nodes that are within a distance of 1. The X-axis shows the number of it- 
erations, and the Y-axis shows the value of 6. Area colors in the heatmap 
depict the average of the following procedure: randomly select a vector v 
and a (5-bounded sequence around v, run the simulation for the randomly 
selected values and return the L^o distance between the last output vector 
and u (calculated as u = W~^v). The heatmap uses a log log scale. Both 
graphs clearly show that as 6 decreases and the number of iterations in- 
creases, the output of SS-Iterative converges to be bounded by a small 
hypercube around u. 

Note that the unit disc weighted topology matrix is characterized 
by Halloo ~ 0-02, ||-B||oo = 0.97 while the circle graph is characterized 
by Halloo ~ 0-33, 1 = 0.66. As expected, using unit disc topology 
requires a larger number of iterations for convergence (depends on 
In addition, in the unit disc topology the value of 6 has a lesser effect 
on the convergence, due to the value of ||^||oo) which affects the minimal 
radius around the output. Since 1 1^| is smaller in the unit disc topology, 
increasing 6 does not significantly affect the convergence. 
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Fig. 2. Sim. of a Circle graph. Fig. 3. Sim. of a Unit-Disc graph. 



6 Extension to the asynchronous model 



Our second novel contribution is in extending our model to support asyn- 
chronous communications. In a large sensor network, it is unreasonable to 
assume that the sensors operate in synchronous rounds. Furthermore, as 
known from the linear iterative algorithms literature, algorithms usually 
converge in less asynchronous rounds (when compared to synchronous 
rounds). 

When considering the asynchronous model, it is more convenient to 
discuss shared-memory as means of communication.^ Thus, assume that 
for each directed edge between Pi,pj there is a read- write register Rp^^p^ 
that is written by pi and read by pj . 

An asynchronous run is an infinite sequence of configurations Cq — >• 
Ci ... such that some process p performs an atomic step between 
configuration Ci and Cj+i. An atomic step consists of reading or writing 
from a single register. Notice that in the current model a configuration 
consists of all of the registers and of the local variables at the different 
nodes. 

In this section we again prove that starting from an arbitrary config- 
uration, when the inputs are bounded, the outputs are bounded as well. 
We consider each configuration Cr to be assigned a vector input I{r) such 
that if node pi reads the input when performing an atomic step on Cr it 
reads the value of Ip^{r). Equivalently, the output vector of configuration 
Cr is 0{r). 



^ In [7] it is shown how to convert an algorithm based on shared-memory to a message- 
passing algorithm with links of bounded capacity. 
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Figure 4 outlines Async-SS-Iterative which is a direct translation 
of SS-Iterative to the shared-memory model. 



Algorithm Async-SS-Iterative 

01 : Forever do: /* executed on node pi */ 

/* write current value of Op- to all neighbors */ 
02: for each pj e iV(p,) 
03: write Op. to Rpi,p^\ 

/* update Op. according to values of neighbors * / 
04: set Op. ;= Wp-^p. ■ Ip.; 
05: for each pj £ N{pi): 
06 : read Rp . ,p . into temp; 

07: update Op- := Op. + Wp.^p^ ■ temp; 

08: od. 

Fig. 4. A self-stabilizing iterative algorithm for asynchronous networks. 



Async-SS-Iterative consists of two phases: in the first, the previous 
value of Op. is written to all its neighbors. In the second phase pi calculates 
its new value of Op. by reading the registers of all its neighbors. 

We consider only "fair" runs, in which each node performs an atomic 
step infinitely many times. Thus, each node performs both phases in- 
finitely many times. A round is defined to be the shortest prefix of a run 
such that each node has performed lines 02-07 in the algorithm. We num- 
ber each atomic step and each round. Note that a round consists of many 
atomic steps. 

We model a fair run as follows. Each node pi performs infinitely many 
atomic steps, and participates in infinitely many rounds. Notice that the 
registers pi reads in round k + I have all been last written to, no earlier 
than during round k. Since a round consists of each node performing 
all the steps in the algorithm, each node pi manages to read all of its 
neighboring registers and write to all of its neighboring registers every 
round. Thus, there is some atomic step r (during round k+l) such that: 

Op^(r) = Wp^^p^ ■ IpXr') + ^ Wp^^p. ■ Op. (r'j) , 

where r' , r'- (for all pj ^ pi) are smaller than r and are from at least round 
k. 
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Let u be such that u = Av + Bu, and let the inputs be from a 5- 
bounded input sequence around v. Denote c(r) = 0(r) — u and z = 
maxj |cp.(0)|. 

Theorem 2. Given a normalized diagonally dominant and invertible W, 
and while considering only fair runs, there are constants ci , C2 , where ci > 
0, and 1 > C2 > 0, such that Async-SS-Iterative e-always converges 
with e{At,6,C) = 6 ■ ci + (02)^* • z; where At counts the asynchronous 
rounds of a fair run. 

Proof. Notice that if pi did not perform the rth atomic step then Op- (r) = 
Op^{r — 1) and therefore Cp. (r) = Cp. (r — 1). Consider the value of Cp. (r) 
when Pi did perform the rth atomic step (during round A: + 1). 

CpAr) = Op, (r) - Up. 

= Wpi,pi ■ Ip, {r') + Wp^^p^ ■ Op^ (r'j) - Wp^^p^ ■ Vp^ - ^ Wp^^p^ ■ Up, 

= Wpi,P. ■ ilpAr') - VpJ + ^ w;p„p, • (Op. {r'j) - Up J 
= Wpi,pi ■ (Ipiir') - VpJ + 5^ Wp,,p,. • Cp.{r'j) , 

where r' and the different r'- are smaller than r and are all from round k 
or round k + \. 

By using Lemma 1 we get: 

|cp,(r)| < |«;p^,p, • (/p,(r') - VpJ| + max |cp^. (r^-)l |w^p.,p, I 

Pj ... 

< ■ {Ip^{r') - VpJI + \\B\\^ \Cp^aA'^max)\ 

< 5 + \^Pmax{'''max')\ 1 

for some Pmax and r^ax < that is from round k ox k + 1. 

Therefore, for any pi during round k + \ there is a list of length 
I > k oi nodes pi,p2, ■ ■ ■ ,P£ and a sequence of length £ of atomic steps 
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ri > r2 > ■ ■ ■ > Vi = 0, such that 



Cp^{r)\<6+\\B 
<S+\\B 
= S-il + 




< S 



Eii^ii- + ii^ii-i^p.(^^)i 



z=0 



5 




Denote by ci = 



1 



, and C2 = Halloo- We have that for node pi 



l-\\B\\ 



OO 



performing the rth atomic step during round k it holds that |cp. (r)| < 



In fair runs, there are infinitely many rounds k, thus, as / and r go to 
infinity, we have that ||0(r)||;^ is bounded by a hypercube of length 5 ■ c\ 
around u. 

7 Discussion 

We have shown that the algorithm SS-Iterative is a modification of the 
Jacobi iterative method to solve a set of equations = b, where A is 
given and b is dynamically changing but bounded. Moreover, Theorem 1 
is a generalization of previous analysis of Jacobi's convergence. Our moti- 
vation for SS-Iterative originates from the sensor calibration problem 
where sensors need to calibrate their noisy readings. Unlike previous ap- 
proaches to this problem, we assume a dynamic system with an infinite 
execution of the algorithm. In this setting the readings of the sensors con- 
tinuously change. Under the assumption that the readings' changes are 
bounded, we have shown that the calibrated output is bounded as well. 

Further application for SS-Iterative can be found in any setting 
where it is desired to solve Tlx = b in a converging and self-stabilizing 
manner, while A is given, and b may change slightly from one round 
to the next. Notice that the analysis given in Section 4 holds in such a 
system. 

As noted in Remark 2 the matrix A is "part of the code" . An optional 
alternative to the current solution is to compute A~^ (the inverted matrix 
of A) beforehand and include it "as part of the code". Thus, each node 
could locally solve x = A~^b, and it can be shown that x will be bounded 



b ■ c\ + c^ - z <b ■ S^- z. 



□ 
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(as long as b is bounded). The main problem with such a solution is the 
connectivity requirements it incurs. In our solution, scalar values are sent 
in the network only between direct neighbors. The matrix W represents 
a weighted adjacency graph. Once inverted, the matrix might not 

be sparse. A non-zero entry w^j^ € means that node pi needs to 

communicate with node pj. This extra communication might cause the 
algorithm to lose its self-stabilizing properties, as non-neighboring nodes 
would require a self-stabilizing overlay network for their communication. 

The assumption of a predefined A is suitable for static networks in 
which the communication graph is predetermined. For dynamic networks, 
it would be interesting to adjust SS-Iterative to discover the connectiv- 
ity of the network, inferring the optimal weights dynamically. We assume 
that after the weights are calculated, the topology of the sensor network 
remains stable, thus the convergence analysis of Section 4 should hold. 

7.1 Relation to Perturbation Theory 

A large amount of research focused on the problem of solving = b 
when A and b are not exactly known. That is, let ^4 = j4 + 6A and 
b = b + (5b, and consider the equation Ax = b; what can be said about 
X in relation to x? 

Our setting is "easier" in one sense, and "harder" in a different sense. 
In our setting A is known, i.e., 6 A = 0. However, b is not well defined. 
That is, the input vector - which is described by b - changes over time. 
When solving Ax = b it is assumed that there is some b that is constant 
but it was measured with an error. In our case, b is not constant as it 
changes over time, while it represents the measurements correctly. 

As a future research, it would be interesting to consider the implica- 
tions of adding inaccuracy to the measurements. The vast body of knowl- 
edge regarding perturbation theory would definitely aid in this extension 
to our model. 

7.2 Relation to convex optimization 

Many practical optimization problems are given in the quadratic form 
f{x) = l/2xAx — b-^x, where the task is to compute miux /(x) distribu- 
tively over a communication network. A survey showing several appli- 
cations can be found in [3]. Example applications are monitoring, dis- 
tributed computation of trust and ranking of nodes and data items. 

A standard way for solving miux /(x) is by computing the derivative 
and comparing it to zero to get the global optimum. When the matrix A 
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is symmetric, /'(x) = j4x — b = 0, and we get a linear system of equations 
A:>c = b. In other words, the convex optimization problem is reduced into 
a solution of a linear system of equations. 

Interior point methods [4, Ch. 11] solve linear programming problems 
by applying Newton method iteratively. Each computation of the Newton 
step involves a solution of a linear systems of equations. An area of future 
work is to examine the applications of our self-stabilizing algorithm to 
these methods. The difficulties arise from the fact that the matrix A needs 
to be recomputed between iterations, so nodes need to be synchronized 
and aware of the current iteration taking place. 
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